logo of company

Intersections of Adversity & Neurodiversity: Adverse Childhood Experiences’ Association to Mental Health & the Buffering Role of Flourishing

This study extends analyses conducted by Bethell et al. (2023), who explored the role of family resilience and connection in promoting flourishing among U.S. children facing adversity. Our study applies a similar framework but focuses on neurodivergent children, examining how adverse childhood experiences (ACEs) influence mental health outcomes and how flourishing behaviors—such as curiosity, emotional control, and task persistence—may buffer these effects. This study was conducted by the Slopen Laboratory at the Harvard T.H. Chan School of Public Health.

Author: Adrián Medina

Date: October 5, 2024

Project Overview

This repository contains files relevant to a study investigating the relationship between adverse childhood experiences (ACEs) and mental health outcomes in children with neurodevelopmental disorders (NDDs), such as Autism Spectrum Disorder (ASD) and ADHD. The study specifically examines the moderating role of child flourishing behaviors—curiosity in learning, emotional control, and task persistence—on the effects of ACEs on mental health outcomes.

Study Design:
The study utilizes data from the National Survey of Children’s Health (NSCH), focusing on children aged 6-17 with reported neurodevelopmental issues (N = 44,776, MAge = 12.2). This cross-sectional analysis examines ACE exposure and its impact on three primary mental health outcomes: anxiety, depression, and behavioral issues. The moderating role of child flourishing behaviors is assessed using interaction terms in logistic regression models.

Data Collection:
Key variables collected include:

  • Adverse Childhood Experiences (ACEs): Derived from parent-reported experiences such as exposure to violence, household dysfunction, and discrimination.
  • Mental Health Outcomes: Parent-reported diagnoses of anxiety, depression, and behavioral problems, confirmed by health care providers or educators.
  • Child Flourishing Metrics: Based on parent-reported behaviors such as curiosity in learning, emotional control, and task persistence, which were used to develop a Child Flourishing Index (CFI).

Analytic Approach:
The primary analytic approach involved logistic regression models to assess the association between ACEs and mental health outcomes, with interaction terms to explore the moderating effect of child flourishing. Key covariates include age, sex, race/ethnicity, and socioeconomic status.

Goals:
The study aims to:

  • Examine the dose-response relationship between the number of ACEs and risks of mental health challenges in neurodivergent children.
  • Identify the protective role of child flourishing behaviors, offering insights into resilience-promoting interventions for this vulnerable population.
Tip

For detailed information on the dataset variables, refer to the NSCH Data Dictionary.

Code Workflow

  1. Data Frame Initialization
  2. Analytic Data Preparation
  3. Frequencies & Descriptive Statistics
  4. Inferential Statistics

Data Frame Initialization

Set up the R environment by configuring the CRAN repository, installing essential packages, setting the base path, and loading the primary dataset into the local environment.

Expand Code 
# Set CRAN repository for consistent package installation
options(repos = c(CRAN = "http://cran.rstudio.com/"))

# Install and load the pacman package for efficient package management
if (!require(pacman)) install.packages("pacman")
library(pacman)

# Use p_load function to install (if necessary) and load packages
p_load(dplyr, tidyr, tidyverse, knitr, ggplot2, Hmisc, broom, stats, kableExtra, webshot2, sjPlot)

# Specify the 'base_path' where you can find your data files, ASSUMING they're in the same directory, & set it as WD
base_path <- "~/GitHub/Adversity-Neurodiversity/data_files"
setwd(base_path)

# Load primary data file
neurodivergent_data <- read_csv("neurodivergent_data.csv")
1
See details under the important callout below!
Data Extraction & Filtering (Archived Code)

The file being used for these analyses is a subset of a “master” file containing the compiled contents of the data releases from the National Survey on Child Health’s 2016-2022 cycles. The master file is nearly 500 MB (0.5 GB) in size, making it cumbersome to process in real-time. For efficiency and ease of analysis, I am using a smaller subset that contains only the relevant data. This subset was generated using the data extraction and filtering process detailed in the archived code below, in an effort to maintain transparency and ensure reproducibility of the workflow.

Expand Code 
# Define the path to the dataset
data_path <- "~/Downloads/Data2_2016to2022.csv"  
Data2_2016to2022 <- read.csv(data_path)

# Prepare project data: Select specific variables, create 'neurodivergent' variable, & filter data based on age and neurodivergent status
neurodivergent_data <- Data2_2016to2022 %>%
  select(year, fpl, fpl_mean, sc_age_years, sc_hispanic_r, sc_race_r, sc_sex, higrade_tvis, k2q35a, k2q35c, k2q38a,
         k2q38c, k2q30a, k2q30c, k2q36a, k2q36c, k2q60a, k2q60c, k2q37a, k2q37c, downsyn, downsyn_desc, 
         k2q31a, k2q31c, k2q61a, cerpals_desc, k2q33a, k2q33b, k2q33c, k2q32a, k2q32b, k2q32c, k2q34a, k2q34b, k2q34c, 
         ace1, ace3, ace4, ace5, ace6, ace7, ace8, ace9, ace10, ace11, ace12, k6q71_r, k7q85_r, k7q84_r) %>%
  # Creating 'neurodivergent' variable as NSCH does not explicitly inquire about neurodivergent status.
  # Using a list of diagnoses provided by NSCH to define neurodivergent individuals.
  mutate(neurodivergent = if_else(k2q35a == 1 | k2q38a == 1 | k2q36a == 1 | k2q60a == 1 | 
                                  k2q37a == 1 | k2q30a == 1 | downsyn == 1 | 
                                  k2q31a == 1 | k2q61a == 1, 1, 0)) %>%
  # Filter data to include only individuals aged 6 or above and identified as neurodivergent
  filter(sc_age_years >= 6, neurodivergent == 1)

# Resulting filtered data to be used for further analysis
write.csv(neurodivergent_data, "neurodivergent_data.csv")

Analytic Data Preparation

Set up recoding of variables by transforming ‘ace1’ into a dichotomous variable, calculating the Childhood Flourishing Index (CFI), and categorizing total ACE counts. Recode education, income, sex, race/ethnicity, and CFI into categorical factors for analysis.

Expand Code 
# Recode 'ace1' to dichotomous variable based on analytic specifications
neurodivergent_data$ace1_recode <- ifelse(neurodivergent_data$ace1 %in% c(2, 3, 4), 1, 
                                          ifelse(neurodivergent_data$ace1 == 1, 0, NA))

# Recode ACE counts and calculate Childhood Flourishing Index (CFI)
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Total ACE count for each child excluding missing values
    ace_total = rowSums(select(., ace1_recode, ace3:ace12) == 1, na.rm = TRUE),
    
    # Categorize total ACEs for further analysis
    ace_counts = cut(ace_total, breaks = c(-1, 0, 1, 3, Inf), labels = c('0', '1', '2-3', '4+'), right = TRUE),
    ace_counts = factor(ace_counts, levels = c("0", "1", "2-3", "4+")), # Recreate 'ace_counts' factor

    # Calculate the Childhood Flourishing Index (CFI)
    CFI = rowSums(select(., k6q71_r, k7q85_r, k7q84_r) == 1, na.rm = TRUE),
    
    # Recode CFI into dichotomous variable
    CFI_dich = case_when(
      CFI == 3 ~ "Flourishing",
      CFI %in% 0:2 ~ "Not Flourishing",
      TRUE ~ NA_character_
    ),
    CFI_dich = factor(CFI_dich, levels = c("Not Flourishing", "Flourishing")) # Recreate 'CFI_dich' factor
  )

# Recode binary outcomes for anxiety, depression, and behavioral issues into factors
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Recode Anxiety (k2q33b)
    Anxiety = factor(k2q33b, levels = c(2, 1), labels = c("No", "Yes")),
    
    # Recode Depression (k2q32b)
    Depression = factor(k2q32b, levels = c(2, 1), labels = c("No", "Yes")),
    
    # Recode Behavioral Issues (k2q34b)
    Behavioral_Issues = factor(k2q34b, levels = c(2, 1), labels = c("No", "Yes"))
  )

# Recode education and income categories
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Recode education into categorical factors
    highgrade_tvis_cat = case_when(
      higrade_tvis == 1 ~ "Less than high school",
      higrade_tvis == 2 ~ "High school (including vocational)",
      higrade_tvis == 3 ~ "Some college or associate degree",
      higrade_tvis == 4 ~ "College degree or higher",
      TRUE ~ NA_character_
    ),
    highgrade_tvis_cat = factor(highgrade_tvis_cat, levels = c("Less than high school", "High school (including vocational)", "Some college or associate degree", "College degree or higher")),

    # Recode federal poverty level categories into factors
    fpl_cat = case_when(
      fpl %in% c(50:99) | fpl_mean < 100 ~ "Less than 100%",
      fpl %in% c(100:199) | fpl_mean < 200 ~ "100% to 199%",
      fpl %in% c(200:299) | fpl_mean < 300 ~ "200% to 299%",
      fpl %in% c(300:399) | fpl_mean < 400 ~ "300% to 399%",
      fpl %in% c(400:999) | fpl_mean < 999 ~ "400% or greater",
      TRUE ~ NA_character_
    ),
    fpl_cat = factor(fpl_cat, levels = c("Less than 100%", "100% to 199%", "200% to 299%", "300% to 399%", "400% or greater"))
  )

# Recode sex and race categories
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Recode sex into categorical factors
    sc_sex_cat = case_when(
      sc_sex == 1 ~ "Male",
      sc_sex == 2 ~ "Female",
      TRUE ~ NA_character_
    ),
    sc_sex_cat = factor(sc_sex_cat, levels = c("Male", "Female")),

    # Recode race and ethnicity
    race = case_when(
      sc_race_r == 1 ~ "White",
      sc_race_r == 2 ~ "Black",
      sc_race_r == 3 ~ "American Indian or Alaska Native",
      sc_race_r %in% 4:5 ~ "Asian, Native Hawaiian, or Pacific Islander",
      sc_race_r %in% 6:7 ~ "Other or mixed race",
      TRUE ~ NA_character_
    ),

    ethnicity = case_when(
      sc_hispanic_r == 1 ~ "Hispanic/Latino",
      sc_hispanic_r == 2 ~ "Non-Hispanic/Latino",
      TRUE ~ NA_character_
    ),

    # Recreate 'sc_race_cat' from race and ethnicity combinations
    sc_race_cat = case_when(
      race %in% c("White") & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic White",
      race %in% c("Black") & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic Black or African American",
      race %in% c("American Indian or Alaska Native") & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic American Indian or Alaska Native",
      race %in% c("Asian, Native Hawaiian or Pacific Islander") & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic Asian, Native Hawaiian or Pacific Islander",
      race %in% c("Other or mixed race") & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic Other or mixed race",
      ethnicity == "Hispanic/Latino" ~ "Hispanic or Latino of any race",
      TRUE ~ NA_character_
    ),
    sc_race_cat = factor(sc_race_cat, levels = c("Non-Hispanic White", "Non-Hispanic Black or African American", 
                                                 "Non-Hispanic American Indian or Alaska Native", 
                                                 "Non-Hispanic Asian, Native Hawaiian or Pacific Islander", 
                                                 "Non-Hispanic Other or mixed race", 
                                                 "Hispanic or Latino of any race"))
  )
1
ace1 is originally coded by NSCH as a frequency measure, but it’s needed as indicative variable of presence like the other ace# variables.
2
Part of this includes the calculating of the Childhood Flourishing Index (CFI) based on Bethell et al. (2019) criteria

Frequencies & Descriptive Statistics

Expand Code 
# Creating frequency tables for specific diagnostic outcomes to assess their distribution in the dataset
freq_k2q33b <- table(neurodivergent_data$k2q33b == 1, useNA="ifany")
freq_k2q32b <- table(neurodivergent_data$k2q32b == 1, useNA="ifany")
freq_k2q34b <- table(neurodivergent_data$k2q34b == 1, useNA="ifany")

# Combine frequency tables into a single data frame for easy viewing
freq_table <- data.frame(
  Diagnosis = c("Anxiety", "Depression", "Behavioral Problems"),
  Count = c(freq_k2q33b["TRUE"], freq_k2q32b["TRUE"], freq_k2q34b["TRUE"])
)

# Use knitr::kable to format the frequency table for better readability in the output document
knitr::kable(freq_table, caption = "Frequency of Specific Diagnostic Outcomes")
Frequency of Specific Diagnostic Outcomes
Diagnosis Count
Anxiety 13725
Depression 6289
Behavioral Problems 13535
Expand Code 
# Count the number of subjects with a CFI score of 3 using a logical condition, as this is considered 'Flourishing.'
num_subjects_cfi_3 <- sum(neurodivergent_data$CFI == 3, na.rm = TRUE)
cat("Number of Subjects Flourishing (CFI Score of 3):", num_subjects_cfi_3, "\n")
Number of Subjects Flourishing (CFI Score of 3): 3265
Expand Code 
# Histogram of ages for neurodivergent participants
ggplot(neurodivergent_data, aes(x = sc_age_years)) +
  geom_histogram(binwidth = 1, fill = "blue", color = "black") +
  labs(title = "Age Distribution of Neurodivergent Participants", x = "Age (Years)", y = "Count") +
  theme_minimal()

Expand Code 
# Bar plot of ACE counts categories
ggplot(neurodivergent_data, aes(x = factor(ace_counts), fill = factor(ace_counts))) +
  geom_bar() +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  labs(x = "ACE Counts Category", y = "Number of Children", fill = "ACE Category",
       title = "Distribution of Children by ACE Counts") +
  scale_x_discrete(labels = c('1' = "0 ACEs", '2' = "1 ACE", '3' = "2-3 ACEs", '4' = "4+ ACEs")) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Expand Code 
# Bar plot of CFI counts categories
ggplot(neurodivergent_data, aes(x = factor(CFI), fill = factor(CFI))) +
  geom_bar() +
  geom_text(stat = 'count', aes(label = ..count..), vjust = -0.5, color = "black") +
  labs(x = "Childhood Flourishing Index (CFI) Score",
       y = "Number of Children",
       fill = "CFI Score",
       title = "Distribution of Children by CFI Score") +
  scale_x_discrete(labels = c('0' = "0", '1' = "1", '2' = "2", '3' = "3")) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Expand Code 
# Bar plot of CFI_dich counts categories
ggplot(neurodivergent_data, aes(x = CFI_dich)) +
  geom_bar(aes(fill = CFI_dich)) +
  geom_text(stat = 'count', aes(label = ..count..), vjust = -0.5) +
  scale_fill_manual(values = c("0" = "blue", "1" = "green")) +
  labs(title = "Distribution of CFI Dichotomous Variable",
       x = "Child Flourishing Index (Dichotomous)",
       y = "Count") +
  theme_minimal() +
  theme(legend.title = element_blank()) # Remove the legend title if not needed

Expand Code 
# Bar plot for educational attainment counts categories
ggplot(neurodivergent_data, aes(x = highgrade_tvis_cat)) +
  geom_bar(aes(fill = highgrade_tvis_cat)) +
  geom_text(stat = 'count', aes(label = ..count..), vjust = -0.5) +
  scale_fill_brewer(palette = "Pastel1") +
  labs(title = "Distribution of Adults' Highest Educational Attainment",
       x = "Educational Attainment",
       y = "Count") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.title = element_blank())

Expand Code 
# Bar plot for sex counts distribution
ggplot(neurodivergent_data, aes(x = sc_sex_cat)) +
  geom_bar(aes(fill = sc_sex_cat)) +
  geom_text(stat = 'count', aes(label = ..count..), vjust = -0.5) +
  scale_fill_manual(values = c("Male" = "lightblue", "Female" = "pink")) +
  labs(title = "Distribution by Sex",
       x = "Sex",
       y = "Count") +
  theme_minimal() +
  theme(legend.title = element_blank())

Expand Code 
#Bar plot for SES counts categories
ggplot(neurodivergent_data, aes(x = fpl_cat)) +
  geom_bar(aes(fill = fpl_cat)) +
  geom_text(stat = 'count', aes(label = ..count..), vjust = -0.5) +
  scale_fill_brewer(palette = "Pastel2") +
  labs(title = "Distribution of Socioeconomic Status",
       x = "Socioeconomic Status",
       y = "Count") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.title = element_blank())

Expand Code 
# Bar plot for race counts categories
ggplot(neurodivergent_data, aes(x = sc_race_cat)) +
  geom_bar(aes(fill = sc_race_cat)) +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5)) +
  scale_fill_brewer(palette = "Set3") +
  labs(title = "Distribution by Race",
       x = "Race",
       y = "Count") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.title = element_blank())

Expand Code 
# Calculate mean and standard deviation of age for neurodivergent participants
age_stats <- neurodivergent_data %>%
  filter(neurodivergent == 1) %>%
  summarise(
    Mean_Age = mean(sc_age_years, na.rm = TRUE),
    SD_Age = sd(sc_age_years, na.rm = TRUE)
  )

# Print the statistics using a nicely formatted table
knitr::kable(age_stats, caption = "Age Statistics of Neurodivergent Participants")

# Calculate prevalence (%) of diagnoses among subject sample
variables <- c("k2q31a", "k2q35a", "k2q38a", "k2q30a", "k2q36a", "k2q60a", "k2q37a", "downsyn", "k2q61a")
variable_names <- c("ADHD", "Autism/ASD", "Tourette’s Syndrome", "Learning Disability", 
                    "Development Delay", "Intellectual Disability", 
                    "Speech/Other Language Disorder", "Down Syndrome", "Cerebral Palsy")

# Pre-calculate total respondents to avoid repeated computation
total_respondents <- colSums(!is.na(neurodivergent_data[variables]))

# Vectorized prevalence calculation
prevalences <- colSums(neurodivergent_data[variables] == 1, na.rm = TRUE) / total_respondents * 100
prevalence_table <- data.frame(Diagnosis = variable_names, Prevalence = prevalences)

# Display prevalence data in a table
knitr::kable(prevalence_table, caption = "Prevalence of Neurodivergent Conditions")

# Prepare the data for cross-tabulations
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    k2q32b = factor(k2q32b, levels = c(2, 1), labels = c("No", "Yes")),
    k2q33b = factor(k2q33b, levels = c(2, 1), labels = c("No", "Yes")),
    k2q34b = factor(k2q34b, levels = c(2, 1), labels = c("No", "Yes")),
    CFI_dich = factor(CFI_dich, levels = c(0, 1), labels = c("Not flourishing", "Flourishing")),
    ace_counts = factor(ace_counts, levels = c(1, 2, 3, 4), labels = c("0", "1", "2 or 3", "4+"))
  )

# Create the cross-tabulations
table_data <- neurodivergent_data %>%
  select(k2q32b, k2q33b, k2q34b, CFI_dich, ace_counts) %>%
  pivot_longer(cols = c(k2q32b, k2q33b, k2q34b, CFI_dich), names_to = "Condition", values_to = "Status") %>%
  group_by(Condition, Status, ace_counts) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(Condition, ace_counts) %>%
  mutate(percent = n / sum(n) * 100) %>%
  pivot_wider(names_from = ace_counts, values_from = c(n, percent), names_sep = ", ")

# Format the table for presentation using 'kable' and 'kableExtra'
knitr::kable(table_data, format = "latex", caption = "Frequency and Percentage of Psychopathology and Flourishing across ACEs Categories") %>%
  kableExtra::kable_styling("striped", full_width = FALSE)

Inferential Statistics

Expand Code 
# Logistic regression models without interactions
model_anx <- glm(Anxiety ~ ace_counts + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years,
               data = neurodivergent_data, family = binomial)

model_dep <- glm(Depression ~ ace_counts + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years,
               data = neurodivergent_data, family = binomial)

model_beh <- glm(Behavioral_Issues ~ ace_counts + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years,
               data = neurodivergent_data, family = binomial)

# Display the models using tab_model with 95% confidence intervals
tab_model(model_anx, model_dep, model_beh, show.ci = 0.95)
1
The tab_model function from the sjPlot package automatically exponentiates the logistic regression coefficients into Odds Ratios when displaying models with a binomial family (logistic regression). 
Profiled confidence intervals may take longer time to compute.
  Use `ci_method="wald"` for faster computation of CIs.
Profiled confidence intervals may take longer time to compute.
  Use `ci_method="wald"` for faster computation of CIs.
Expand Code 
# Plot for the Anxiety model with x-axis limits between 0.1 and 5
plot_model(model_anx, type = "est", show.values = TRUE, show.p = TRUE, ci.lvl = 0.95, 
           title = "Odds Ratios - Anxiety", axis.lim = c(0.1, 5))
Profiled confidence intervals may take longer time to compute.
  Use `ci_method="wald"` for faster computation of CIs.
Expand Code 
# Plot for the Depression model with x-axis limits between 0.1 and 5
plot_model(model_dep, type = "est", show.values = TRUE, show.p = TRUE, ci.lvl = 0.95, 
           title = "Odds Ratios - Depression", axis.lim = c(0.1, 5))
# Plot for the Behavioral Issues model with x-axis limits between 0.1 and 5
plot_model(model_beh, type = "est", show.values = TRUE, show.p = TRUE, ci.lvl = 0.95, 
           title = "Odds Ratios - Behavioral Issues", axis.lim = c(0.1, 5))
Profiled confidence intervals may take longer time to compute.
  Use `ci_method="wald"` for faster computation of CIs.

  Anxiety Depression Behavioral Issues
Predictors Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p
(Intercept) 11.02 6.82 – 18.33 <0.001 6.93 4.08 – 11.99 <0.001 49.42 34.07 – 72.76 <0.001
ace counts [1] 1.29 1.09 – 1.53 0.003 1.24 1.02 – 1.51 0.033 1.14 0.99 – 1.31 0.063
ace_counts2-3 1.52 1.28 – 1.81 <0.001 1.47 1.22 – 1.78 <0.001 1.12 0.98 – 1.29 0.095
ace counts [4+] 1.88 1.54 – 2.29 <0.001 1.80 1.47 – 2.20 <0.001 1.36 1.17 – 1.59 <0.001
sc race cat [Non-Hispanic
Black or African
American]		 0.76 0.59 – 1.01 0.053 0.94 0.72 – 1.25 0.658 1.03 0.86 – 1.24 0.742
sc race cat [Non-Hispanic
American Indian or Alaska
Native]		 0.52 0.31 – 0.93 0.018 1.35 0.68 – 3.06 0.430 0.91 0.57 – 1.52 0.707
sc race cat [Non-Hispanic
Other or mixed race]		 0.89 0.71 – 1.14 0.349 1.01 0.80 – 1.29 0.917 0.93 0.78 – 1.10 0.387
sc race cat [Hispanic or
Latino of any race]		 0.70 0.59 – 0.84 <0.001 0.84 0.69 – 1.02 0.071 0.90 0.78 – 1.05 0.169
sc sex cat [Female] 1.35 1.19 – 1.53 <0.001 1.42 1.25 – 1.61 <0.001 1.21 1.09 – 1.35 0.001
fpl cat [100% to 199%] 0.95 0.76 – 1.19 0.670 0.88 0.70 – 1.10 0.273 0.91 0.77 – 1.09 0.310
fpl cat [200% to 299%] 1.03 0.81 – 1.30 0.799 0.82 0.65 – 1.04 0.109 0.83 0.69 – 0.99 0.039
fpl cat [300% to 399%] 0.84 0.66 – 1.07 0.152 0.78 0.61 – 1.01 0.057 0.69 0.57 – 0.83 <0.001
fpl cat [400% or greater] 0.91 0.72 – 1.15 0.443 0.68 0.54 – 0.86 0.002 0.71 0.59 – 0.84 <0.001
highgrade tvis cat [High
school (including
vocational)]		 1.12 0.72 – 1.69 0.605 1.10 0.72 – 1.65 0.651 1.00 0.72 – 1.36 0.981
highgrade tvis cat [Some
college or associate
degree]		 1.11 0.72 – 1.65 0.627 1.14 0.74 – 1.68 0.539 0.98 0.71 – 1.33 0.887
highgrade tvis cat
[College degree or
higher]		 1.22 0.79 – 1.83 0.353 1.04 0.68 – 1.54 0.864 0.86 0.62 – 1.16 0.330
sc age years 0.97 0.95 – 0.99 0.001 0.96 0.93 – 0.98 <0.001 0.85 0.84 – 0.87 <0.001
Observations 14612 7372 15339
R2 Tjur 0.007 0.018 0.040
Expand Code 
# Logistic regression model with CFI_dich as a moderator for Anxiety
model_anx_mod <- glm(Anxiety ~ ace_counts * CFI_dich + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years,
                     data = neurodivergent_data, family = binomial)

# Logistic regression model with CFI_dich as a moderator for Depression
model_dep_mod <- glm(Depression ~ ace_counts * CFI_dich + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years,
                     data = neurodivergent_data, family = binomial)

# Logistic regression model with CFI_dich as a moderator for Behavioral Issues
model_beh_mod <- glm(Behavioral_Issues ~ ace_counts * CFI_dich + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years,
                     data = neurodivergent_data, family = binomial)

# Display the models with tab_model
tab_model(model_anx_mod, model_dep_mod, model_beh_mod, show.ci = 0.95, transform = NULL)
1
Since we included interaction terms, we want to see the raw coefficients instead. You can set transform = NULL in the tab_model function to ensure the values aren’t exponetiated into Odds Ratios. 
Profiled confidence intervals may take longer time to compute.
  Use `ci_method="wald"` for faster computation of CIs.
Profiled confidence intervals may take longer time to compute.
  Use `ci_method="wald"` for faster computation of CIs.
Expand Code 
# Plot marginal effects for the Anxiety model
plot_model(model_anx_mod, type = "int", title = "Interaction: ace_counts * CFI_dich (Anxiety)")
# Plot marginal effects for the Depression model
plot_model(model_dep_mod, type = "int", title = "Interaction: ace_counts * CFI_dich (Depression)")
# Plot marginal effects for the Behavioral Issues model
plot_model(model_beh_mod, type = "int", title = "Interaction: ace_counts * CFI_dich (Behavioral Issues)")

  Anxiety Depression Behavioral Issues
Predictors Log-Odds CI p Log-Odds CI p Log-Odds CI p
(Intercept) 2.45 1.97 – 2.96 <0.001 1.96 1.43 – 2.51 <0.001 3.95 3.58 – 4.34 <0.001
ace counts [1] 0.26 0.08 – 0.43 0.004 0.20 -0.00 – 0.40 0.053 0.14 -0.00 – 0.28 0.053
ace_counts2-3 0.41 0.23 – 0.58 <0.001 0.34 0.15 – 0.53 0.001 0.12 -0.02 – 0.26 0.090
ace counts [4+] 0.59 0.39 – 0.79 <0.001 0.54 0.34 – 0.75 <0.001 0.31 0.15 – 0.46 <0.001
CFI dich [Flourishing] -0.96 -1.43 – -0.46 <0.001 -1.43 -2.23 – -0.66 <0.001 -0.76 -1.31 – -0.17 0.009
sc race cat [Non-Hispanic
Black or African
American]		 -0.28 -0.54 – 0.00 0.047 -0.06 -0.33 – 0.22 0.674 0.03 -0.15 – 0.22 0.735
sc race cat [Non-Hispanic
American Indian or Alaska
Native]		 -0.66 -1.18 – -0.07 0.018 0.28 -0.40 – 1.11 0.452 -0.11 -0.57 – 0.41 0.671
sc race cat [Non-Hispanic
Other or mixed race]		 -0.13 -0.36 – 0.11 0.287 0.01 -0.23 – 0.25 0.953 -0.07 -0.24 – 0.10 0.412
sc race cat [Hispanic or
Latino of any race]		 -0.35 -0.53 – -0.16 <0.001 -0.18 -0.37 – 0.02 0.076 -0.10 -0.25 – 0.04 0.165
sc sex cat [Female] 0.31 0.18 – 0.43 <0.001 0.35 0.23 – 0.48 <0.001 0.19 0.08 – 0.30 0.001
fpl cat [100% to 199%] -0.05 -0.27 – 0.18 0.692 -0.12 -0.35 – 0.10 0.296 -0.11 -0.29 – 0.06 0.211
fpl cat [200% to 299%] 0.03 -0.21 – 0.27 0.792 -0.19 -0.43 – 0.04 0.112 -0.21 -0.39 – -0.03 0.024
fpl cat [300% to 399%] -0.17 -0.41 – 0.07 0.174 -0.25 -0.51 – 0.00 0.055 -0.40 -0.59 – -0.21 <0.001
fpl cat [400% or greater] -0.08 -0.32 – 0.15 0.487 -0.38 -0.61 – -0.14 0.002 -0.37 -0.55 – -0.19 <0.001
highgrade tvis cat [High
school (including
vocational)]		 0.08 -0.36 – 0.50 0.700 0.10 -0.33 – 0.50 0.650 -0.02 -0.36 – 0.29 0.880
highgrade tvis cat [Some
college or associate
degree]		 0.06 -0.37 – 0.47 0.762 0.12 -0.30 – 0.52 0.554 -0.05 -0.38 – 0.26 0.754
highgrade tvis cat
[College degree or
higher]		 0.15 -0.29 – 0.56 0.480 0.03 -0.40 – 0.43 0.893 -0.18 -0.51 – 0.12 0.256
sc age years -0.03 -0.05 – -0.01 0.003 -0.04 -0.07 – -0.02 <0.001 -0.16 -0.17 – -0.14 <0.001
ace counts [1] × CFI dich
[Flourishing]		 -0.26 -0.95 – 0.43 0.456 -0.08 -1.21 – 1.07 0.897 -0.47 -1.25 – 0.29 0.227
ace_counts2-3:CFI_dichFlourishing 0.05 -0.64 – 0.75 0.895 1.37 0.30 – 2.51 0.014 -0.50 -1.25 – 0.24 0.188
ace counts [4+] × CFI
dich [Flourishing]		 0.51 -0.48 – 1.67 0.344 1.07 -0.03 – 2.26 0.064 -1.08 -2.05 – -0.13 0.026
Observations 14612 7372 15339
R2 Tjur 0.011 0.023 0.046

This document was made by Adrián Medina as part of the Slopen Laboratory, based on the Lumo quarto template.